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INTRODUCTION 


the stochastic difFeretial equation of the Langevin type 


Gaussian random walks prove to be a natural and 
rather universal starting point for many stochastic pro¬ 
cesses. In fact, the famous central-limit theorem shows 
that many independent random movements of finite vari¬ 
ance tr^ = {x^) always pile up to display a Gaussian dis¬ 
tribution [T]. In particular, Gaussian random walks con¬ 
stitute the basis of the most important tool in the theory 
of financial markets, the Black-Scholes option price the¬ 
ory [3 (Nobel Prize 1997), by which a portfolio of assets 
is hoped to remain steadily growing through hedging [3] . 

However, since the last stock market crash and the still 
ongoing financial crisis it has become clear that distribu¬ 
tions which describe realistically the behaviour of finan¬ 
cial markets belong to a more general universality class, 
the so-called Levy stable distribution They result 

from a sum of random movements of infinite variance [8] , 
and account for the fact that rare events, the so-called 
Black-Swan Events [5], which initiate crashes, are much 
more frequent than in Gaussian distributions. These are 
events in the so-called Levy tails oc l/\x\^^^ of the dis¬ 
tributions, whose description is based on a generalized 
Hamiltonian [10]: 

iL(p) = const (p2)V2, (1) 


^x{s) = x(s) = r](s), (2) 

as 

where r]{s) is a noise variable as a function of a pseu¬ 
dotime s with zero expectation value and a probability 
distribution characterized by a parameter A [18] : 


P[t]] = 


f dsH{rj) _ 


= / Ppexp 


ds 


iprj - 



Using this we may solve the stochastic differential equa¬ 
tion in which the noise 77(5) has nonzero correlation 
functions for even n = 2,4,6 ,... : 


(? 7 (si)... 77(s„)) = J X>?7 77 (si)... 77(s„)P[?7]. (4) 

For A = 2, the distribution is Gaussian and r]{s) is a stan¬ 
dard white noise variable. If we solve ([^ in I? dimensions 
with an initial condition x(0) = 0 , the variable x(s) has 
a distribution 

Pg(x, s) = (47 rs)“'°'^^e“’‘(5) 

This distribution is the Green function of the Fokker- 
Planck equation 


Such tail-events are present in many physical situations, 
e.g., in velocity distributions of many body systems with 
long-range forces HD, in the self-similar distribution of 
matter in the universe [HHH], and in the distributions 
of windgusts m and earthquakes uni, with often catas¬ 
trophic consequences. 

Distributions with Levy tails are a consequence of 
rather general maximal entropy assumptions m- In the 
limit A —> 2, the Levy distributions reduce to Gaussian 
distributions. 

The simplest Levy-type random walk is described by 


(5s+P^)H g(x,s) = (5(s)(5(^^(x), (6) 

where p = iSx = iV. For X ^ 2, the distribution is 
non-Gaussian and it solves the fractional Fokker-Planck 
equation 

[5s + (p2)^/^]P(x,s)=<5(s)5(^)(x). (7) 

A solution of this equation that evolves from the 15- 
function is 


P(x,s) = e-"(p')"^'<5(^)(x) 


( 8 ) 
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and for s = 1 it coincides with the noise probability, 

P(x,l)|..^ = P(r,) = I (9) 

Applications of the fractional Fokker-Planck equation are 
numerous in non-Brownian diffusion processes. These are 
observed in chaotic systems and in the fluid dynamics of 
rheology and biology. See [HI |20] for an overview. The 
mathematics of Eq. 0 with variable diffusion coefficient 
is in pT] . 

The fractional Fokker-Planck equation 0 can be gen¬ 
eralized further to the double fractional Fokker-Planck 
equation 


can be derived from a random walk model when the mean 
waiting time of the walker diverges [28] . 

It is the purpose of this note to calculate the Green 
functions of general fractional Fokker-Planck equation 
(101 and specify the path integrals solved by them [HEl]. 


DOUBLE FRACTIONAL FOKKER-PLANCK 
EQUATION 

A convenient definition of the fractional derivatives 
uses the same formula as in the dimensional continua¬ 
tion of Feynman diagrams [301 El] > 

(p2)A/2 ^ j (14) 


where p 4 = dt, p = id^ = *V and a parameter has been 
allowed for that is the analogue of the diffusion constant 
D in the ordinary diffusion process [22] . 

We should explain the physical origin of the fractional 
powers in the space and time derivatives of the above 
equation. Such powers occur naturally in many-particle 
systems if the interaction strength or the range becomes 
very large. As long as the interaction strength is small 
and the range is short, such systems are described by a 
second-quantized field theory with a free-particle action 

Ao = Jdtd^xtp^ (x, t){idt-\-h'^\7'^/2m — V{x))ip{x, t), (11) 
and an interaction of the type 

Aint = ^ J dtd^x {'if if . (12) 

The partition function can be calculated from the func¬ 
tional integral 


The solution of ( |10[ ) can be written formally as 

P(x,t) = [(p4 + e)i-^+PA(p2)"/2]-'<5(t)d(^)(x), (15) 

where inhnitesimal e > 0 ensures forward-in-time nature 
of the Green function, and its explicit appearance will 
be suppressed from now on. Using the representation 
6{t) = we arrive at 


= J 


dE 


^ — iEt 


2tT (-iP)l-7+DA(p2)V2 


SP^x). (16) 


Now we expand the fraction into a geometric series, and 
integrate term by term using the formula [32] 


^+oo 


dE 


„-iEt 


e{t)t' 


i(l-7)-7 


J_^ 2tt {-iE + e)iP-A(n+i) r[(l-7)(n-hl)]’ 

(17) 

where 9{t) is the Heaviside step function. The result can 
be cast as 


Z = j, VifVif'<e^^-^°+-^'''^P^. (13) 

A perturbation expansion leads to an effective action in 
the form of a power series of where 4' = {if) are the 

expectation values of the held. This series is divergent 
and must be resummed. For large interaction strength 
(/, this produces anomalous power behaviors in the held 
strength as well as in the momenta [23l|24|- The free-held 
part of the effective action leads to a held equation of the 
fractional Fokker-Planck or Schrddinger type, in which 
momentum and energy appear with powers diherent from 
A = 2 and 7 = 0, respectively. 


In addition, equations of the type (10) are known to 


govern various different phenomena. In chaotic systems, 
for example, they describe anomalous diffusion processes 
with memory (time non-locality) [23 [20|- In fact, the 
fractional time derivatives also arise as the inhnitesimal 
generators of coarse grained time evolutions or they 


P{x,t)=e{t)t 

(18) 

where Ea,p{z) = T{an+ii) Mittag-Leffler 

function [33] |34| . This can be interpreted by writing 


P(x,t) = (x| ^^.(t) |0), (19) 

with the 7 -deformed evolution dehned by 

U^{t) = 6»(t)t-^Ai_.^,i_^(-A-^i7), (20) 

with H = Ua(p 2)^A gSj (See Fig. 0 ) . The occur¬ 
rence of the Mittag-Leffler function in solutions of the 
time-fractional Fokker-Planck equation has been noted 
previously, for example, in the review article [22] . 


For 7 = 0, the equation (10) reduces to a single (space) 


fractional Fokker-Planck equation 

\p4 + Dx{p‘^)^^'^]P{x,t) = S^^^{x)S{t), (21) 











Pt{0-5,s) 
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Figure 1: (Color online) The function U-y{t) for H — 1, and 
various values of 7 . Dotted (blue) curve: 7 = 0, standard 
exponential function; Dashed (red) curve: 7 = 0.1; Solid (yel¬ 
low) curve: 7 = 0.5. 



Figure 2: (Color online) Dotted (blue) curve: A = 2, stan¬ 
dard Gaussian distribution; Dashed (red) curve: A = 1.5; 
Solid (yellow) curve: A = 1, Cauchy-Lorenz distribution. The 
length scale is is = 2{DxsYP. 


the Mittag-Leffler function reduces to Ei^i{z) = exp( 2 ;), 
and the evolution operator recovers its standard form 
?7o(^) = 0{t) exTp{—tH). The solution, which we shall 
denote by Px (x, t) for a more specific reference, is the 
multivariate Levy stable distribution [36) : 


Px{x,t) 


-»px 

(27r)^ 


( 22 ) 


For A = 2, it reduces to the standard quantum mechani¬ 
cal Gaussian expression (§. For A = 1, the result is 


[r(i :>/2 -b l/2)/7r(^+i)/2]£)^i 

[{Dxty + \x\^]D/2+1/2 ’ 


(23) 


which is the Cauchy-Lorentz distribution function. In 
Fig.@ we plot Px in Z? = 1 dimension for A = 1,1.5, 2. 

In the Appendix we provide various useful representa¬ 
tions of Px(x, t). At this place it is worth mentioning 
that this probability can be written as a superposition 
of Gaussian distributions PG{a,x.) = {4:TTa)~^^‘^e~^ 
to be specified in Eq. ([^. 



Ft(1, s) 


-Pt(1.5,s) 


s 

2.0 


(a) 


(b) 


Figure 3: (Color online) (a) PT{t,s) as a distribution of t 
with increasing values of the pseudotime s = 0.5,1,1.5. (b) 
PritjS) as a distribution of s with increasing values of the 
real time t — 0.5,1,1.5. In both cases 7 = 0.03. 


Smeared-time representation, and relation between 
physical time t and pseudotime s 


lods 


If we use in (16) the Schwinger’s formula l/A = 


, we can express P(x, f) as an integral 

POO 

P{x,t) = / dsPx{x,s)PT{t,s), 

Jo 


(24) 


where Px solves the space-fractional diffusion equation 


( 21 ), with t = s, and Pt solves the time-fractional equa¬ 


tion 

-\-pl~'^]PTit,s) = S{t)S{s), (25) 

which encodes the relation between the pseudotime s 


and the physical time t. The factorized ansatz (24) has 


been used previously in m to solve the time-fractional 
Fokker-Planck equation. 


For 7 = 0 , PT{t,s) = S{t — s), and (24) reduces to 
P(x,t) = Pxix,t). 

For 7 > 0, we obtain an asymmetric Levy stable dis¬ 
tribution |38] 


PTit,s) = [ 
J —( 


dE 

2t: 




(26) 


An important feature is that Prit, s) vanishes for t < 0. 
This can be seen by placing the branch cut of a multi¬ 
valued function along the negative real axis, and 

calculating ( [ 2 ^ as a complex integral with contour that 
follows the real axis, and closes in the upper half-plane. 
See Fig. (a) where Pt is plotted as a function of t for 
the case 7 = 0.03, and various values of s. 


It is illustrative to view formula (24) as a smearing 


of the distribution Pjf(x, s) around the time position t, 
defined by the probability density function PT(t, s). For 
this purpose we plot in Fig. (b) Pt(Z, s) as a function 
of s, with parameter t describing the position of the peak 
in the probability distribution. 

The two plots in Fig. ([^ are related through the for¬ 
mula 


PT{t,s) = {C/t)PT{C,C^-^E-h), 


(27) 
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Figure 4: (Color online) PT{t,s) as a function of both t and 
s. Here 7 = 0.1. 


which can be deduced from (26) by a simple change of 
the integration variable E —)■ {C/t)E. Here C is an ar¬ 
bitrary constant. The function PT{t,s) as a function of 
two variables is shown in Figure (|^. 

When 7 = 0, Frit, s) = 6{t — s) is concentrated at the 
point t, i.e., there is no smearing. For increasing 7 the 
peak around t broadens, which can be accounted for by 
derivatives of the 5-function. The action of Pt on a test 
function /(s) is 


PT{t,s)f{s) = 

We represent 
calculate 



“ roo 

- - ds PT{t,s){s-t)^. 

u=0 -^0 

(28) 

= (-1)”/ dT5^'^'>{T - t)f{T), and 



x|/^t 


Figure 5: (Color online) In all cases A = 2. Dotted (blue) 
curve: 7 = 0, standard Gaussian distribution; Dashed (red) 
curve: 7 = 0.03; Solid (yellow) curve: 7 = 0.1. 


Fox H-function representation of the Green function 


Solution of the double fractional equation (10 1 has 
been obtained previously in terms of the Fox H-function 
|39) . We derive the same result starting from for¬ 
mula (24), where we consider the representation (58) of 
Px(x, s). Integration over the pseudotime s can be per¬ 
formed, followed by the E integration, that yields 


P(x,t) =- 


f-T' 


.t3/2|xP 2.3 


( 

■|xr 

X 


jt. 



(l,l);(l-7,l-7) 


(l.l),(i3/2.A/2);(l.A/2)y 


(33) 


Here £* = 2{Dxt^ 7 ) 1 /a jg g, t-dependent length scale, 
and iF-function nasi], defined by the 

contour integral 


n T, u ^ k fdE fc!6»(t)t(i-'>')'=-T' 

j^dsPT[t,s)s -J — 


to find that 


ffl 

PT{t,s) = —^Cn{t)S^^\s - t), 

ni 


n—O 


where 


■w-i: H'-i) 


k=0 


k\9{t)t 


--i(k+l) 


r[(i- 7 )(fc + i)]- 


(29) 

(30) 

(31) 


dz r(i + z)r(f + ^z)T{-z) 


F^(x,t)|xp 

t-iTr-D/2 7^ 27rf r(-|z)r(l-7-b (l-7)z)L (-i 

(34) 

where the contour C runs from —zoo to -|-zoo. In Fig. 

we show how values of 7 > 0 modify the Gaussian 
distribution (for which A = 2, 7 = 0). 


The large-|x| asymptotics of (33) is governed by the 
pole of the integrand at 2 : = 1 : 


P|x|^P(x,t) 


|x|->.oo 


-r(^) 


7rW2r(2-27)r(-f) 


(35) 


In view of these relations, the equation (24) translates 
into 

°° 1—fi" 

^(x,t) = y'— l-Cn{t)d^Px{y^,t). (32) 

n\ 

n=0 


Analysis of the small-|x| behavior is more subtle due 
to a richer pole structure of the integrand in (34) (see 
[42]). If we assume only simple poles, we can extract the 
leading behavior 


and 


Ixl-i-n 


One can easily verify that for 7 = 0, c„ = 5 „ 0 j 
P(x,t) = Px(x,t). 


PP(x, t) 


A(t)-fP(t)|x|2A-^, 2X-D<2 
A(t)-f C>[|x|2](t), 2X-D>2 ^ ’ 
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Figure 6: (Color online) Dotted (blue) curve: 7 = 0, A = 
1, Cauchy-Lorentz distribution; Dashed (red) curve: 7 = 
0.1, A = 1; Solid (yellow) curve: 7 = —0.1, A = 1, assumes 
negative values. Here D=3. 


with 

A'= J drlipx- (43) 

Now P{x,t) = (xt |0 0 ) satisfies the ordinary fractional 
Fokker-Planck equation 

[P4 + Dx{p^)’'/^]P{x,t) = (44) 

which has been discussed at length in recent literature 

mi- 


LANGEVIN EQUATIONS AND COMPUTER 
SIMULATIONS 


with 


A{t) = 


B{t) = - 


^i-D/ 2 ^-C 2 /y 

sin(7r'f)ixf)r[ (^-^)f-^) ]’ 

7r-f^/2r(f - A) 


r(A)r(7-i)^2A- 


(37) 


(38) 


In particular, for 2X < D the value of P{x,t) tends to 
either +00 or —00 as |x| —^ 0. See Fig. 


PATH-INTEGRAL FORMULATION 


We note that the probability (15) may be calculated 


from the doubly fractional canonical path integral over 
fluctuating orbits t(s),x(s) p 4 (s),p(s) viewed as func¬ 
tions of some pseudotime s [i5] : 


{xbtbSb\xatsSa} = J VxDt'Dp'Dpie-^, 


(39) 


with A being the euclidean action of the paths t(s),x(s): 

.4 = y ds[i{px'— ip 4 t') —'H{p,P 4 )]. (40) 

Here t'{s) = dt{s)/ds, x'(s) = dx{s)/ds, and 'H{p,p 4 ) = 
p\~"' + D\{p^)^^'^. At each s, the integrals over the com¬ 
ponents of p(s) run from —00 to 00 , whereas those over 
P 4 {s) run from —zoo to zoo. To obtain the distribution 
P{x,t), we finally form the integral 


poo 

P{x,t) = / ds{xts| 000 }. 

Jo 


(41) 


This is analogous to prescription (24) which links solu¬ 


tions of the space- and time-fractional diffusion equations 


(21) and (25). 


If 7 = 0, the path integral over P 4 {s) yields the func¬ 
tional 5[t'(s) — 1 ], which ensures that dt and ds increments 
are equal. This brings (391 to the canonical path integral 


{Xbtblxja) = j VxDpe-^ , 


(42) 


In the past, many nontrivial Schrodinger equations 
(for instance that of the 1 /r-potential) have been solved 
with path integral methods by re-formulating them on 
the pseudotime axis s, that is related to the time t via 
a space-dependent differential equation t'{s) = f{x{t)). 
This method was invented by Duru and Kleinert |45j to 
solve the path integral of the hydrogen atom, and has re¬ 
cently been applied successfully to various Fokker-Planck 
equations gi El- The stochastic differential equation 
(47), that connects pseudotime s and the physical time t, 


may be seen as a stochastic version of the Duru-Kleinert 
transformation that promises to be a useful tool to study 
non-Markovian systems. 


Certainly, the solutions of Eq. (44) can also be obtained 


from a stochastic differential equation 


X = ?7, 


(45) 


whose noise is distributed with a fractional probability 


PH = 




(46) 


Simulating this stochastic differential equation on a com¬ 


puter, we confirm the analytic form ( 22 ) of Px(x, s) = 
P(x, t) for 7 = 0 . See Fig. ^ (a). 


Analogously, the solution of Eq. (25) can also be ob¬ 
tained from a SDE 


i'{s) = t]t{s), 
with noise distribution 


P[qT,] = J P, 


(47) 


(48) 


and compared with the result (26) for PT{t,s). See Fig. 

[!l(b). 

Solution of the double fractional Fokker-Planck equa- 
pa can be obtained, in view of the relation (|4T|) 


tion 


(or (24)), by simulating (45) for t = s and (47), and let¬ 


ting the final value of the pseudotime s be random. This 
yields a probability distribution P{x,t). In Fig. 1^ we 


























6 



Figure 7: (Color online) Comparison of analytic (solid red 
curve) and numerical (blue circles) results for the distribution 
function Px{x, s=l) in D = 1 dimension (a), and PT{t,s=l) 
for 7 = 0.3 (b). In each case an average has been taken 
over 5000 representative trajectories of stochastic differential 
equations (|45|l and (471, with 10 time steps As = 0.1. 



Figure 8: (Color online) Comparison of computer simulation 
and the renormalized exact solution P(x, t) for t=0.2,l. 


compare the results of a computer simulation with the 
analytic form (34) by plotting P(x, t) as a function of x 
for various values of time t. Since the distribution P(x, t) 
itself is not normalized, but rather 


J d^xP{^, t) = J^dsPrit, s) = (49) 

we define a renormalized version P(x, t) = 
P(x, t)// d^xP{x,t). 


SUMMARY 

Summarizing, we have seen that a many-body system 
with strong couplings between the constituents satisfies a 
more general form of the Schrodinger equation, in which 
the momentum and the energy appear with a power dif¬ 
ferent from A = 2 and 7 = 0, respectively. We have 
calculated the associated Green functions and discussed 
their properties and their representations. We pointed 
out that these Green functions can be written as path 
integrals over fluctuating time and space orbits that are 


functions of some pseudotime s. This is a Markovian 
object, but non-Markovian in the physical time t. The 
non-Markovian character is caused by the fact that func¬ 
tion t{s) follows a stochastic differential equation of the 
Langevin type. 

The particle distributions can also be obtained by solv¬ 
ing a Langevin type of equation in which the noise has 
correlation functions whose probability distribution is 
specified by an equation like (46). 

The Green functions whose theory was presented here 
will play an important role in the development of an in¬ 
teracting theory of fields whose worldlines contain non- 
Gaussian random walks displaying extremely large devi¬ 
ations from their avarages. 


Acknowledgment: We are grateful to P. Jizba, and 
A. Pelster for useful comments. V.Z. received support 
from GACR Grant No. P402/12/J077. 


Appendix 1: Fractional differential operators that 
enter the general fractional Fokker-Planck equation ([^ 
are defined through formula (14). Using (5(x) = 
{4:TTa)~^^'^e~^ and = 5{t — a), we derive 

the following relations, 


Ixl^ = 






= r(a-bl)(p4) “ S{t), 


(51) 


which we can substitute into ( |33t, ( |34[ ) in order to verify 
that these satisfy the equation (10). We first obtain 


4^(x, t) = J^ ^r(i + z)r(-z)ui(p2)^-/2 

X(p4)(7-1)G+I)j(0)(x)^(i), 

which can be pole-expanded to yield 


(52) 


n—0 

Summing up this geometric series, we arrive at 

P(x,t) = [pI-^ + (54) 

Appendix 2: We derive several expressions for the so¬ 
lution Px(x, s) of (21), starting from the representation 

On expanding the exponential, and representing the 
powers as (p2)W2 = r[-An/2]-i 
the momentum integration yields the superposition of 
Gaussian expression 


Px{yi,s)= ( —fx{a)PG{y^,D^^^s'^/^a), 

Jo ^ 


( 55 ) 
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with weight 


= Y. 


h "!r(-An/2) ■ 


(56) 


To prove this, we perform the cr-integration term by term, 
using the formula /“ = T{i/)/a^, and ob¬ 

tain the large-|x| expansion 


Px(x,s) = 


rD/2| 


\D 


E 


(_l)nr(An^) 


tiZ'o ^!r(-An/2) 


(57) 


where ig = . The series can also be viewed as 

a pole expansion of the contour integral, and hence 


1 

f dz V{^)V{-z) 

[1x^1 

TT-^/^jx]-® J 

Q 27rz r(—Az/ 2 ) 

cc 


Px(x,s) = 

7r^/“|X|^ fr i\—/\z/z; ^,v 

\58) 

with the contour C running from —ioo to -|-*oo. From 
this, the expansion ( [57| arises by enclosing the right com¬ 
plex half-plane and calculating the residua of the inte¬ 
grand, using Res(r(az -I- 6), —(n -I- h)/a) = (—l)"'/(n!a). 
A small-|x| expansion of (58) is obtained by closing the 
integration contour in the left half-plane, leading to 


P lx . v (-i)-2/A r(gi^) 

^ > 2^ ^D/ 21 D „!r(f+n) 


n=0 


121 


(.1 


(59) 


The series (57) and (59) are convergent, or asymptotic, 
or even trivially zero, depending on the parameter A. 
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